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Abstract 

The products of the two and three bipolar harmonics 3^£j^f (1*31 , ^2) are represented as the finite 
sums of powers of the three relative coordinates r^^r^i and r%\- The complete (angular+radial) 
integrals of the products of the two and three bipolar harmonics in the basis of exponential ra- 
dial functions are expressed as finite sums of the auxiliary three-particle integrals r n & /(a, /3, 7). 
The formulas derived in this study can be used to accelerate highly accurate computations of the 
rotationally excited (bound) states in arbitrary three-body systems. In particular, we have con- 
structed compact variational wave functions for the triplet and singlet 2P(L = 1)— states in light 
two-electron atoms and ions. 
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The bipolar harmonics y^M( r 3U r 32) [it are extensively used in various methods devel- 
oped for highly accurate solutions of different three-body problems originated in atomic, 
molecular and nuclear physics. The functions, Eq.(j2J), are often used to represent the 'an- 
gular dependence' of the wave functions. A very high efficiency of the bipolar harmonics 
as 'angular' functions follows from their explicit form which reflects a number of simple 
and transparent physical ideas. There are a number of successfull generalizations of bipolar 
harmonics to the four- and five-body systems, where the three- and four-polar harmonics 
arise. In this study, however, we restrict our analysis to the three-body systems only. In 
general, the three-body bipolar harmonics are written in the form 

3^M 2 (x,y)= E ^w^iW^iWy) (i) 

mi 7712 

where 

\i2m2 are Clebsh-Gordan coefficients, while x and y are the two vectors 
uniformly related with the three-body system. Let us designate three particles in such a 
system by the numbers 1, 2 and 3. The 'natural' choice of the two fundamental vectors 
for an arbitrary three-body system is x = r-3 — ri = and y = r3 — r 2 = r 32 . In this 
case each bipolar harmonic, Eq.([T]), takes the form 3^lm ( r 3i> r 32) an d contains the vectors 
r 31 and r 32 . Here and everywhere below in this study we assume that these two vectors are 
truly independent, i.e. r 31 ^ Ar 32 , where A is a numerical constant. In actual cases there 
are some additional advantages, if one operates with the angular functions of the unit norm 
vectors n 3 i and n 3 2, where ny = ^ M -. In these vectors the bipolar harmonics are written in 
the from 

3^|(r 31 ,r 32 ) = 4\ril £ C^ i£2m Y hmi (n 31 )Y hm2 (n 32 ) = 4r|^(n 31 , n 32 ) (2) 

mim,2 

Variational wave functions which include bipolar harmonics are used to aproximate the 
actual wave functions for the bound states with the non-zero angular momentum L. In 
many papers the bipolar harmonics are called and considered as the 'angular parts' of basis 
functions. It is assumed that the additional 'radial' part of the total wave function depends 
upon the three radial coordinates r 32 ,r 31 and r 21 only, i.e., it does not contain any of the 
angular variables. It is clear that the bipolar harmonics with the same L and M values (or 
indexes) form the (2L + 1)— dimensional representation of the rotation group 50(3). The 
explicit form of the matrixes which describe transformations of the bipolar harmonics during 
rotations can be found with the use of Eq. (j2J) and formulas from Q and 3. It can be shown 
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that each matrix element of these matrix is the product of the two Clebsh-Gordan coefficients 
and two Wigner's D— functions. The explicit expression for these matrix elements can be 
reduced to another 'short' form, but below we will not need these formulas. 

In this communication we develop the new method to operate with the bipolar harmonics. 
Our main interest below is related with the products of the two and three bipolar harmonics 
and angular integrals of such products. First, note that bipolar harmonics with the same LM 
indexes form the closed algebra, i.e. the product of two bipolar harmonics 3^lm ( r 3i> r 32) an d 
y l LM ( r 3i) r 32) is always represented as the finite sum of bipolar harmonics with the same 
values of L (angular momentum of the system) and M (its z— projection). This can be 
written in the following form 

^(r 31 ,r3 2 )^(r 31 ,r3 2 ) = £ f LM {li,h, 4, 4; 4, W£$(r 3 i, r 32 ) (3) 

where /£m(^i, ^2> ^3> 4> 4) -4) are the numerical coefficients which can be determined, e.g., 
by integration of the both parts of Eq. over all angular variables. 

Let us brielfy discuss the integration over the angular variables in an arbitrary non- 
relativistic system of three particles. In general, to describe the non-relativistic quantum 
system of spinless particles one needs 9 (3 x 3) dynamical variables. Three of these nine 
variables describe the translations of the solid triangle of particles. The internal state of 
three-body system does not change during such translations. In actual cases these three 
(Galilean) translations can be separated by using an appropriate choice of internal coordi- 
nates. Formally, we can assume that after such a separation of translations one of the three 
particles, e.g. the third particle, will always be at rest. The remaining six coordinates are 
separated into two groups: (a) three coordinates which are rotationally invariant, i.e. they 
do not change during any rotation of the whole three-body systems, and (b) three coor- 
dinates which describe rotations of the whole three-body system. The scalar coordinates 
which do not change during any rotation and/or translation of the three-body system can 
be chosen as the three interparticle distances r 32 ,r 3 i,r 2 i. 

The choice of the three truly independent 'rotational' coordinates in the three-body sys- 
tem is slightly more complicated, since such coordinates must be related with the angular 
coordinates of the two vectors r 3 i and r 32 . Let 8i, <pi and # 2 , <p2 be the shperical coordinates 
of these two vectors. In these coordinates for an elementary volume dV we can write 

dV = r^dr^sind^^rl^dr^sinB^m (4) 
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These six coordinates (r 3 i, #31, 3i , r 32 , 6 32 , ^32) can be used to describe an arbitrary three- 
body system. However, as it was shown by Hylleraas in 1929 4J it is better to chose 
three-body coordinates in a different way. In [4] three radial variables were chosen as scalar 
interparticle distances (or interparticle coordinates) =| rj — |= r^, where (ij) = (32), 
(31), (21). On the other hand, it is clear that three angular variables can be chosen as the 
three Euler's angles 03i, #31, 032. In reality, one finds a number of advantages in calculations, 
if such a 'natural' choice of angular variables is used. In particular, in all earlier papers 
(see, e.g., jsj], j(| and references therein) three radial r 3 2,r 3 i,r2i and three angular variables 
03i) $3i) 032 were chosen as described here. In these variables an elementary volume dV takes 
the form 

dV = r 32 r 31 r 21 dr 32 dr 31 dr 21 sin9 31 d9 31 d(f) 32 d(f) 31 (5) 

This part of our study can be concluded with the two following comments. First, the 
three radial r 32 , r 3 ±, r 2 i and three angular variables 03i, #31, 032 are semi-separated from each 
other. This means that that angular integral of any function of these six variables inte- 
grated over three angular variables (or over Euler's angles) is the function of the three radial 
variables only. Furhtermore, in many actual cases these functions of radial variables are 
written in a relatively simple, finite-term form. Second, the three radial variables r 32 , r 31 , r 21 
are not independent of each other, since e.g., r 2 i < r 32 + r 3 \ and r 3 \ >| r 32 — r 2 i |. Such 
constraints substantially complicate analytical and numerical calculations of three-body in- 
tegrals. Therefore, in actual cases it is better to use three truly independent perimetric 
coordinates Ui,u 2 ,u 3 Q], which are related with the relative coordinates by the following 
linear transformation 

^1 = ^31+^21-^21) , u 2 = -(r 32 + r 21 -r 3 i) , u 3 = ^(r 32 + r 31 - r 21 ) 

The inverse relation takes the form 

^32 = u 2 + u 3 , r 31 = ui + u 3 , r 32 = Ul + u 2 

The Jacobian of the (r 32 ,r 31 ,r 21 ) — > (ui,u 2 ,u 3 ) transformation equals 2. Three perimetric 
coordinates ui,u 2 , u 3 are independent of each other and each of them varies between and 
+00. 



4 



The formula for the angular integral of the product of the two bipolar harmonics can be 
written in the form 



(6) 



The explicit form of the radial F^ i2 . iaii function can be found with the use of Eq.© (see, 
e.g., [5] and references therein). The result is 

Fl x l 2 ;tzlSr?>1, ^31, r 2l) = -(' 



31 



1 h h A i 


'441] 


v o o o) < 


£2 ii A J 



[4][4][4][4]E(-i)^l 

A 



Px(x) 




X 



(7) 



where [a] = 2a + 1 and the notation P\(x) stand for the Legendre polynomial of the order A, 
where A is an integer positive number. Also, in this formula (and in some formulas below) 
we use the standard notations for the 3j— and 6j— symbols |8|. The sum over A in Eq.(J7J) 
is always finite, since the product of two 3j— symbols is not zero only for those A which are 
bounded btween the following values: max{\ £\ — £ 3 |, | ^2 — ^4 |} < A < min{£i + £ 3 , £ 2 + £4}. 
Moreover, the product of these two 3j— symbols equal zero unless the two sums of the 
corresponding momenta [£\ + £ 3 + A and £2 + £4 + A) are even numbers. 
The variable x, in Eq.flTJ), is the following dimensionless ratio 

2 I /v>2 ,v, 2 

x = T ^ + ^ 2 !J1 /gx 

2r 31 r 32 V ; 

This expression can be transformed with the use of the formula (8.911) from [9J for the 
Legendre polynomial P\(x) 



Px{x) 



where A 



1 * (-l) k (2X-2k)\ 
¥^ k\(X-k)\(X-2k)\ 



is the integer part of | and coefficients a\ t u are 



x 



X~2k 



1 



X-2k 



(9) 



k=0 



'-l) k (2X-2k)\ 



k\(X - k)\(X - 2k)\ 

Now, by using the formula, Eq.flB]), one finds the following expression for the x x ~ 2k factor 
from Eq.flHD 

X-2k X-2k n 



X 



X-2k 



„2k-X„2k-X 
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31 
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where k < A (see Eq.([9l)) and notation C% stands for the binomial coefficients (the number 
of combinations from b by a, where a and b are positive integer numbers). 

The formula, Eq. fllQI) . allows one to re- write the expression, Eq.(jH]) in the form 

i A A— 2k n 

TD („\ — \^ ,, \ " r< n ( 1 \-2k-2n 2k-\+2n-2m 2m 

^v^) 1J\ 2—/ \—2k v I Ly n r 32 r 31 r 21 K 11 ) 

Z k=0 n=0 m=0 

Now, by using the formulas, Eqs. (fT0|) and ( TTTj) . we can derive the following finite-sum ex- 
pression for the exponential integral of the Legendre polynomial P\(x) 

Ie 2+tr! £ 1+h (a,l3,r,P\) = J J J p x(x) exp(-ar 32 - /3r 31 - jr 2 i)4l +u+1 rl\ +£3+1 r 21 dr 32 dr 3 idr 2 i 

A A-2fc n 

— 7^\2^2 a ^ k 2^2 ^\-2k 2~L ( — ^) m C n l re 2+ e i+ \-2k-2n+i,e 1 +e :j +2k-\+2n-2m+i,2m+i{a, (3, 7) (12) 

Z k=0 n=0 m=0 I 1 

where I\; in (a, b, c) are the basic three-body integrals defined in jlo|. The definition of the 
basic three-body integral is written in the form 

r fc;Z;n (a,/3,7) = / / / r^r^r^ exp(-ar 32 - f3r 31 - -yr 21 )dr 32 dr 31 dr 21 (13) 



where all indexes k, I, n are assumed to be non-negative integer numbers. The analytical 
formula used for numerical computations of such integrals is obtained from Eq. fl 1 3 1) by 



performing integration in perimetric coordinates js] 

r ( n r ^ = 9 V V V r k r l r n + (k-h + m)] (n- ni + h)\ 

k ' l '' n{ ' P ' 11 k % ^ n % (a + /3y-h+ki+l ( a + T )fc-fci+»i+l (/? + T )n-n 1+ J 1+ l 

k I n (~ik\ s~il\ rm\ 

= 2 A- 1 /' n\ Y^ Y^ Y^ "-"i+^i ^fc-fci+ii ^i-h+m 

' ' '^o6>n^o(« + /?)^ 1+fe1 ^ 

where C]JJ are the binomial coefficients. The formula, Eq. (jl4p . can also be written in a few 
other equivalent forms. The function x ^ +1 in Eq.( TT4"]) is the A n (X) function introduced by 



Larson 11 L For the first time, I have produced the formula, Eq. (I14jl . in the middle of 1980's 
6] and references therein). The formula, Eq. (jl4j) . has been used in calculations of 



(see, e.g. 

various three-body integrals, e.g., integrals containing one or two Bessel functions [10 ] . 

With the use of the formulas derived above one can obtain the closed (i.e. finite term) 
analytical formula for the following exponential integral 

1 
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dr 31 dr 21 = -(-1)VNN[4]NE(-1) A ^ 

A 

A A-2fc n 

2^2 a \k 2^2 ^\-2k zJ ( _ l) m CTr^+^4+A-2A:-2n+l/i+£3+2fc-A+2n.-2m+l,2m+l(a, b, c) (15) 
k=0 n=0 m=0 
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The derivation of this formula was the main goal of our study. This formula is of great 
interest for numerical calculations of matrix elements which are needed to determine the 
total energies of bound states in three-body systems with L > 1 and calculate various 
expectation values. Recently, we have developed a number of fast numerical approaches to 
calculate the auxiliary three-particle integrals r n> ^(a, b, c). This allows one to accelerate 
substantially all numerical calculations of matrix elements (see below). 

Note that the matrix elements of the potential energy are written in the form of Eq. (|14p 
only in those cases when all interparticle interaction potentials are the scalar functions of 
interparticle distances r 32 ,r 31 and r 2 i- Such cases include Coulomb three-body systems and 
three-body systems in which the potential energy is written in the form V a {r^2) + H(^3i) + 
Vc( r 2i)- In more complex cases the interaction potential between each pair of particles can 
also be a function of angular coordinates. Analytical formulas for the matrix elements in 
such cases must include angular integrals from the products of three bipolar harmonics. 
Such integrals are discussed below. 

Now, let us present present the results of numerical calculations with the use of formulas 
derived above. These results are shown in Tables I and II. In Table I we demostrate the 
results of numerical calculations of the l£ 2+ £ 4 .£ 1+ g 3 (a, P, 7; Px) and 

^ 3 £ 4 ( a ) c ) integrals, 

Eq. lfTSl) . determined for different numerical values of its arguments (a,b,c) and parameters 
(L, (^1,^2), (^3,^4))- The method of numerical computations of these (exponential) integrals 
is based on the formula, Eq. (fl5l) . However, for actual calculations this formula has been 
modified to avoid numerical instabilities which arise during summation of large numbers of 
positive and negative terms. The first formula, Eq. (II 21) . was re- written in the form 

A A-2A: n 

2 A 



where A 



k=0 n=0 m=0 

(16) 

, while the coefficients | \ are 

{2\-2k)\ 



k\(X - k)\(X - 2k)\ 

The sums of the positive and negative terms in Eq. (fl6|) must be calculated separately. At 
the second step we have used the following formula 

/ 



£2 U A 




x 



£2 £1 A 

These two formulas are used for very fast and accurate calculations of the exponential 
integrals J r e J 1 e 2 -e 3 e 4 ( a ^ ^ c ) which can be found in each matirx element of the Hamiltonian 
and overlap matrices for the bound states with L > 1. By performing extensive numerical 
computations we have found that the mehtod based on the modified formulas, Eqs. ffTB"]) - 
(JT7]) . is fast, numerically reliable and can be applied in computations of different rotationally 
excited states, including highly excited states with L > 15 — 20. 

In Table II we determine the total energies of the bound P(L = 1)— states of a number 
of three-body systems. These systems include the two-electron °°He, 4 He and 3 He atoms 
and two-electron He-like ions: Li + , Be 2+ , B 3+ and C 4+ . For the bound P(L = 1)— states 



we can perform highly accurate numerical calculations by using our old approach and 
new method described in this study. Therefore, we can compare the final accuracy of both 
methods and computational times needed to compute the same values. The nuclear masses 



of the 3 He and 4 He nuclei are 5495.8852 m e and 7294.2996 m e , respectively 12]. The masses 

of the nuclei in all two-electron ions are assumed to be infinite. For each of these systems 

we determine the total energies E (in atomic units) of the singlet 2 1 P— states and triplet 

2 3 P— states (see Table II). Our trial wave functions contain N = 400 exponential basis 

functions. The explicit form of such functions is 

1 . N 

^32,r 3 i,r 21 ) = -=[1 + (-l) £ Pi 2 ] ■ yiV 2 (r 3 i,r 32 )^C l ea;p(-a l r3 2 - &r 31 - 7< r 21 ) (18) 
V2 i=1 

where Ci are the linear variational coefficients, e = 1 in the case of the tripet states and 

e = 2 (or 0) in the case of the singlet states. The operator P\i is the permutation of the 

two identical particles (electrons 1 and 2) and iV is the total number of terms in the trial 

function. Analogous wave functions with iV = 700 basis functions will later be used as a 

short-term cluster functions in our highly accurate computations of the bound 2 1 P— and 

2 3 P— states in these atomic systems. All such calculations are performed with the use of 



our two-stage optimization strategy [l3j . [14]. 



Note that the matrix elements of the potential energy are written in the form of Eq. ( |T4"1) 
only in those cases when all interparticle interaction potentials are the scalar functions of 
interparticle distances r 32 ,r 31 and r 21 . Such cases include Coulomb three-body systems and 
three-body systems in which the potential energy is written in the form T4(r 32 ) + V;,(r 3 i) + 
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V c (r2i). In more complex cases the interaction potential between each pair of particles 
can also be a function of angular coordinates. In the general case, such potentials can 
be approximated by the sums of bipolar harmonics with the different L and M values. 
Analytical formulas for the matrix elements in such cases must include angular integrals 
from the products of three bipolar harmonics. 

Let us discuss the formulas for the products of three bipolar harmonics 
^L t ,M ( r 3i ) r 3 2), yi±£i h (Tai,*aa) and y%% B (r 31l r 32 ). The general formulas for such products 
can be found in [l|. It is clear that the angular integral from the product of three bipo- 
lar harmonics must be proportional to the Clebsh-Gordan coefficient C L a M a LcMc , or to the 
corresponding 3jm— symbol (see below). Second, as we have mentioned above the (2L + 1) 
bipolar harmonics 3^,jv!( r 3i) r 32) (with the same L, but different M) are the basis vectors 
of the (2L + l)-dimensional representation of the rotation group. Therefore, as it follows 
from Schur's lemma the angular integral of the product of two bipolar harmonics m 
and y^M b is always zero unless the values of L a , M a and Lb, Mb are exactly the same, i.e. 
L a = Lf, and M a = Mb. This expalins the explicit form of the angular integral used in 
Eq.(j6]). For the product of three bipolar harmonics the situation is different and we cannot 
assume a priori that even some of the indexes are equal. In general, the angular integral 
from the product of three bipolar harmonics is written in the form 



/ dnyi% a (v 31 , r 32 )yt%(r u , r 32 )^l(r 31 , r 32 

/ 

M n . Mb M, 



GiiM^hAiA ( r 3 2 ' r 31> r 2l) 



L a Lb L c 

^M a M b M c 

L " U Lc I J2 b Al +ei+e6 ri\ +i ^P x (x) (19) 



where the function G^^^^^r^, r 3 \, r 2 i) depends upon three relative coordinates r 32 , r 3 i 
and r 21 . The explicit formula for this function is obtained from the last equality in Eq. ffTPT) . 
The coefficients b\ in Eq. fTl9l do not depend upon the relative coordinates, but it is a 
functions of all ten values of angular momenta A, L a , Lb, L c , £±,£2, £3, £4, £5, £%■ The formula 
for these coefficients takes the form 

1 / 

h = — (-l) Lb+A \/ [L a ] [L b } [L c ] [£1] [£ 2 ] [£ 3 ] [£4] [£5} [4] • [A] EE 

Ai A2 
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(20) 



In this form Eqs. (ll9l) - (120]) look very similar to Eq.(J6J). Numerical caculations of the 
b\ coefficients with the use of Eg. (120 p is straightforward. For instance, for 4 = 1,4 = 
3,4 = 2,4 = 2,4 — 2,4 — 2 and for L a = 1,-4 = 2,L C = 1 one finds from 
Eq.fl20D 61 = 0.15921224404155089- lO" 1 ,^ = -0.67600638318413508- 10~ 2 and b 5 = 
0.27451117819426586 • 10~ 2 . The coefficients 6a with other values of A equal zero identi- 
cally. 

Finally, as it seen from the formulas, Eqs. (jl9p - (I2"0"j) and Eq. fllip . the calculation of the 
three-body integralas which contain the products of three bipolar harmonics is reduced to the 
computation of some finite sums of the basic (or auxiliary) three-body integrals r^^a, /3, 7), 
Eq. (jl~4l) . The explicit expression for the radial integral of the G^^^. gi £:j ^(^32, ^31, ^21) 
function is 

^MfcfeA = J J J ^+^+ i (i-A+e 3 +e 5 ( r 32i r 31 , r 21 )r 3 2r3ir 21 dr 3 2dr 31 dr2i 

= ^hlta+U+le-A+ta+lsfaPiriPx) > (21) 
A 

where It 2+ £ 4+ t 6; t 1+ e 3+ e B (a, /3, 7; P\) is defined by Eq. (fTB1) . This is another result which is of 
great interest for highly accurate computations of many actual three-body systems. Our test 
calcualtions of the T^i^e^^ e 3 e 5 coefficient perfomed for 4 = 1,4 = 3,4 = 2,4 = 2,4 = 
2,4 = 2, L a = 1, L b = 2, L c = 1 lead to the following result T^. lj2>2 = 1.9488412125971230- 
10 4 . 

We have considered the products of the two and three bipolar harmonics ^j^(r 31 , r 32 ). 
It is shown that angular integrals of such products are represented as the finite sums of 
powers of the three relative coordinates r 32 ,r 31 and r 2 i (or interpartilce distances). The 
six-dimensional (angular + radial) integrals of the products of the two and three bipolar 
harmonics in the basis of exponential radial functions are expressed as finite sums of the 
auxiliary three-particle integrals T nt k t i(cx, fl, 7). The formulas derived in this study can be 
used to accelerate highly accurate computations of the rotationally excited (bound) states 
in arbitrary three-body systems. The methods developed in this study have been used to 
construct very compact (but highly accurate!) variational wave functions of the triplet and 
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singlet 2P(L = 1)— states in light two-electron atoms and ions. 
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TABLE I: Numerical values of the Ig 2+ g 4 .g 1+ g 3 (a,b,c;P\) and 7^^.^^(0,6,0) integrals computed 
for different values of the £±, £2, £3, £4, A and L parameters, a = 1.55,6 = 1.33 and c = 1.07 in all 
cases. 



£1 


£2 


h 


£4 


A 


hih-M^b, c;P x ) 


L 




1 





1 





1 


0.20963264469930568E+00 


1 


0.13008789710134043E+00 


1 


1 


1 


1 


1 


0.11032101479039074E+01 


2 


0.72370337608449528E+00 


1 


1 


2 





1 


0.11491009319449536E+01 


2 


0.69783486804905589E+00 


1 


2 


2 


1 


2 


0.33053241063178740E+02 


3 


0.83514739108659169E+01 


1 


2 


1 


2 


2 


0.34640478424614411E+02 


3 


0.69505739570548458E+01 


3 





1 


2 


2 


0.18295480775566101E+02 


3 


0.14561365717576028E+02 


3 








3 


3 


0.58890134034481354E+02 


3 


0.29445067017240677E+02 



TABLE II: The total energies E of the bound l P(L = 1)- and 3 P(L = 1) -states of some two- 
electron atoms and ions (in atomic units). The total number of basis functions used to construct 
these short-term wave functions is 400. 





EQ-P{L = l)-state) 


E( 3 P(L = l)-state) 


°°He 


-2.12384 308649 747 


-2.13316 419077 725 


4 He 
3 He 


-2.12354 565412 914 
-2.12344 834501 189 


-2.13288 064210 349 
-2.13278 787470 796 


°°Li+ 

oo Be 2+ 
oog3+ 
ooq4+ 


-4.99335 107777 776 
-9.11077 162291 209 
-14.4772 832652 478 
-21.0933 323133 814 


-5.02771 568139 659 
-9.17497 314302 556 
-14.5731 376921 585 
-21.2217 106964 559 
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